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ABSTRACT 

We present a cross-spectra based approach for the analysis of CMB data at large angular scales 
to constrain the reionization optical depth r, the tensor to scalar ratio r and the amplitude of 
the primordial scalar perturbations With respect to the pixel-based approach developed so 
far, using cross-spectra has the unique advantage to eliminate spurious noise bias and to give 
a better handle over residual systematics, allowing to efficiently combine the cosmological 
information encoded in cross-frequency or cross-dataset spectra. We present two solutions to 
deal with the non-Gaussianity of the Q estimator distributions at large angular scales: the 
first one relies on an analytical parametrization of the estimator distribution, while the second 
one is based on modification of the Hamimache&Lewis likelihood approximation at large 
angular scales. The modified HL method (oHL) is powerful and complete. It allows to deal 
with multipole and mode correlations for a combined temperature and polarization analysis. 
We validate our likelihoods on numerous simulations that include the realistic noise levels 
of the WMAP, Planck-LFl and Planck-HFl experiments, demonstrating their validity over a 
broad range of cross-spectra configurations. 

Key words: cosmological parameters - cosmic microwave background - methods: data anal¬ 
ysis - methods: statistical 


1 INTRODUCTION 

One of the main challenges left for the present and future Cosmic 
Microwave Background (CMB) experiments is the high precision 
measurement of the CMB polarization anisotropies at large angu¬ 
lar scales. This signal is extremely interesting because it encodes 
unique informations about the ionization history of the Universe 
and the inflationary epoch and it can be used as an independent and 
complementary probe to the small scale CMB information to con¬ 
strain two important cosmological parameters: the optical depth to 
reionization t and the tensor-to-scalar ratio parameter r which is 
related to the amplitude of the primordial tensor modes. Moreover, 
the large scales CMB signal is useful in breaking parameter de¬ 
generacies, in particular concerning the two parameters: r and the 
amplitude of the primordial scalar density perturbations A , which 
are strongly correlated through the amplitude of the first acoustic 
peak^A^T- = A^e^^’’. 

Current CMB experiments, in particular the ones that, as 
Planck ( [Planck Collaboration I|20f5) , targeted the accurate mea¬ 
surement of the CMB temperature anisotropies, have now reached 
a level of precision and resolution such that they have exploited all 
their statistical power, and are now limited by the systematic ef¬ 


fects related to the instrument design and technology. An unprece¬ 
dented accuracy and care at each step of the data analysis and its 
interpretation is therefore required to access the cosmological in¬ 
formation encoded in the CMB polarization anisotropies at large 
angular scales. In this paper we address this issue focusing on the 
importance of developing statistical methods specific to the analy¬ 
sis of CMB data at large angular scales that allow to minimize the 
impact of residual systematics related to the experimental configu¬ 
ration and design. 

Given that the distribution of the CMB anisotropies is com¬ 
patible with a Gaussian distribution, all the relevant statistical in¬ 
formation is encoded in the two points correlation function of the 
CMB temperature and polarization anisotropies or, equivalently, 
its projection in harmonic space: the angular power spectrum of 
the CMB temperature and polarization fields. This is defined as 
Cf = {atm<Pt>^>)5tc'^ where af„, are the coefficients of the spherical 
harmonic decomposition. The connection between the measured 
CMB data and the theory is done through the CMB likelihood func¬ 
tion = P(d|Cf(a)) that quantifies the match between the CMB 
data d and a given theoretical model parametrized e.g. by a theoret¬ 
ical power spectrum Cf (tr) defined in terms of a set of cosmological 
parameters a. 
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So far the analysis of the CMB anisotropies at large angular 
scales has mostly been based on methods that relies on low resolu¬ 
tion maps in order to compute the exact CMB likelihood function 
in pixel space, ^ = P{A\Ci(a)), with d h M(p) = 2;™ 

This approach is based on the fact that, given that the CMB 
anisotropies are compatible with a gaussian distribution with ran¬ 
dom phases, the follow a multi-variate Gaussian distribution. 
The likelihood function, written in pixel space or, equivantely, in 
terms of the atm coefficients, is gaussian and therefore can be com¬ 
puted exactly (|Gorski et al.||1994| |Slosar et al.||2004| |Page et al.| 
[20071 |Bennett|2013|l. 

The problem of this approach is that, in the case of a real CMB 
experiment, the maps consist in a combination of signal, noise, in¬ 
strumental systematics and must account for the incomplete sky 
coverage necessary to minimize the impact of the foregrounds con¬ 
tamination. In order to achieve the required accuracy at large an¬ 
gular scales, the noise matrix in pixel space must be reconstructed 
with extremely high accuracy to avoid spurious bias on the param¬ 
eters reconstruction. However this accuracy can be extremely hard 
to achieve given the difficulty of the precise characterization not 
only of the noise but also of the residuals systematics related e.g. to 
the instrument, the scanning strategy and the residual foregrounds. 

Alternatively, the likelihood function could be defined in the 
harmonic space as done e.g. in the small scales analysis where 
the data compression from CMB maps to angular power spectra is 
necessary for computational and numerical reasons. However, the 
complication of working in harmonic space at large angular scales 
(low-f multipoles) is related to the fact that the distribution of the 
Cl estimators at low-f is non-Gaussian. In harmonic space the Cg 
consist in the sum of the square of the harmonic coefficients af„, 
and they have a reduced-;^^^ distribution. Therefore the likelihood 
of a theoretical power spectrum as a function of the measured Q 
is non-Gaussian. Contrary to the small-scales analysis, the CMB 
low-f analysis is particularly concerned by this issue given that the 
central limit theorem cannot be invoked. Previous studies, jPerci-j 
jval & Brown|2006||Hamimeche & Lewis|2008^ , developed a CMB 
analysis on large angular scales based on the likelihood definition 
in harmonic space in terms of auto-spectra, that is to say CMB an¬ 
gular power spectra obtained from a given single frequency/dataset 
CMB map. This approach however shared problems similar to the 
pixel based likelihood approach, in particular in terms of the de¬ 
pendency to the noise and of the accurate characterization of the 
systematics effects at the auto-spectra levels. 

In this paper we propose to extend the cross-spectra based ap¬ 
proach for the analysis of the CMB temperature and polarization 
anisotropies to the large angular scales. We provide different so¬ 
lutions to deal with the non-Gaussianity of the cross-spectra es¬ 
timators at large angular scales. Working in harmonic space us¬ 
ing the cross-spectra allows to get rid of noise biases and to min¬ 
imize the residuals systematics effects by exploiting the cross¬ 
correlation between different CMB maps, e.g. cross-frequency and 
cross-datasets. In this sense, the use of cross-spectra allows to ac¬ 
cess the cosmological information encoded in the CMB maps at 
different frequencies and to combine different CMB datasets in a 
more powerful way with respect to the pixel based or auto-spectra 
approach. 

We present a detailed description of the cross-spectra statis¬ 
tics in Sect.j^ In Sect. we describe the Ci estimator that we use 
for the cross-spectra reconstruction and we define the specifications 
used to generate realistic cross-spectra simulations based on pub¬ 
licly available CMB data. Furthermore, in Sect. |3.2| we present the 
formalism to deal with the non-Gaussianity of the cross-spectra Q 


estimators at large angular scales in the case of our realistic simula¬ 
tion settings. Based on this formalism, in Sect. we construct two 
types of cross-spectra based likelihoods: in Sect ]4. Ij we present an 
analytical solution based on the parametrization of the Ci estima¬ 
tor distribution that is useful for the simplest case of a single-field 
analysis where correlations can be neglected. In Sect. |4.2| we then 
define a more general method that allows to easily deal with a joint 
temperature and polarization analysis accounting for both correla¬ 
tions between multipoles and modes (TE, TB, EB). This more gen¬ 
eral method is based on the extension of the lHamimeche & LewisI 
( [2008^ (H&L) approach to the large angular scales analysis and it 
relies on a re-definition of the H&L variable transformation allow¬ 
ing to approximate the CMB likelihood function by a multivari¬ 
ate Gaussian at low multipoles and for cross-spectra. In Sect. we 
present the likelihood results in the case of a single-field analysis, 
describing the validation tests and a comparison of the different 
methods. As the reference single-field we consider the E-modes 
polarization to constrain the optical depth to reionization param¬ 
eter T. The results for the general modified H&L solution (oHL) 
that accounts for the full temperature and polarization analysis in¬ 
cluding all correlations are described in Sec. [^ where we present 
constraints of the t, r and A , parameters. Also, we discuss the opti¬ 
mality of the oHL method with respect to the pixel based likelihood 
solutions. Einally in Sect.jTjwe present our conclusions. We provide 
in the appendix Sect.|^the details of the analytic description of the 
cross-spectra distribution and in Sect.|^we discuss the comparison 
of the auto-spectra and cross-spectra statistics. 


2 CROSS-SPECTRA STATISTICS 


In order to gain some understanding of the underlying statistics, 
we start by presenting the analytical formalism to deal with CMB 
cross-spectra on the full sky, which will be generalized in Sect.j^in 
particular to a cut-sky. We consider the CMB angular power spec¬ 
trum obtained by combining the harmonic coefficients ai„ of two 
different full-sky maps (A, 6), measured with different noise spec¬ 
tra Nf and A®. For a realistic experimental setting the harmonic 
coefficients are also convolved with the beam functions of the two 
maps A and B, and bf. The cross-spectra statistics is defined as: 


^ AxB 

Lf 


2£ 


1 

TT 2 


( 1 ) 


In the Eq. 0 and in the following we make the hypothesis that the 
noise and the residual systematics are not correlated between the 
maps so that the cross-spectra are unbiased estimate of the CMB 
signal. The cross-spectra distribution is given by (we refer to the 
Appendix Sect.j^for the details of this calculation): 




2 (v-i )/2 ^Y{NI2) ^[z{crAcrB)'^l^ 


( 2 ) 


where c = , z = (1 -p^)crA(TB, N = 2t+\is the number 

of modes, is the modified Bessel function of the second kind 
and order v and: 


<Ts=^C^ + N^ 

p= , ^ =• 

+ NfXCf + Nf) 


(3) 
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Figure 1. Examples of the full-sky cross-spectra distributions = 

for N = 5 modes {£ = 2) with ctaitb = 0.1 and varying the degree 
of correlation (cf. Eq. |^). 


Some examples of the shapes for these distributions are shown 
on Fig. [T] where it is interesting to see how the distribution changes 
when varying the degree of correlation between the two maps (p). 
Note that, unlike for auto-spectra, the p^^{c) function can be neg¬ 
ative since the negative exponential decay is compensated by the 
rise of the Bessel function, especially when the noise is important 
(large cr^, o-fi). 

From the characteristic function Eq. jA12| > we can compute 
the cumulant generating function K(t) = In 0(f) and by Taylor- 
expanding it in powers of {it) around zero we obtain the first cu- 
mulants: 

= Cf (4) 

2{Cff + Cf{Nf+N^) + N^N^ 

KiiCe ) = - - - (5) 

- AXB 8(C*)2 -I- 6Cf{Nf + Nf) + 6NfNf 

kACi ) = Cf - ^^(6) 

This generalize the results from |Hamimeche & Lewis] ( |2008[ Ap¬ 
pendix C)) obtained with identical noise. According to the Central 
Limit Theorem, the cumulants above K 2 disappear with the number 
of modes N and the distribution approaches a Gaussian with a vari¬ 
ance given by Eq. 0- Unlike for the auto-spectrum case (Eq. l |B3^ ), 
the estimator ki does not depend on the noise reconstruction. The 
clear advantage of using the cross-spectra is that the estimator is un¬ 
biased whatever knowledge we have of the noise spectra. Also, the 
statistical loss for using cross-spectra with respect to auto-spectra 
is small and minimized if the noise levels of the two maps involved 
are not too different, as shown in details in Sect. Note that in 
general these conclusions hold true also when an incomplete sky 
coverage is considered. 


3 CROSS-SPECTIfA ESTIMATOR 

As in the auto-spectrum case (e.g. |Wandelt et al.|2001^ , the inclu¬ 
sion of some cut on the sky and of anisotropic noise complicates 
the description of a cross-spectrum estimator by correlating modes 
between them (both in £ and m for a non azimuthal mask) and even¬ 
tually distorting the marginal distributions. We then need to rely on 
realistic simulations to take into account the full complexity of the 
problem. 


EE 



multipole 


BB 



multipole 


Figure 2. Polarized power spectra for E modes (top) and B modes {bottom 
for different reionization histories: late (t = 0.056, dashed line) or early 
(t = 0.09, solid line). The primordial B-modes spectra are shown for r = 0.2 
(solid and dashed thick lines) and the tensing contribution to the B-modes 
signal is shown as the thin line. The noise levels for the four considered data 
cases are over plotted {WMAP, Planck-10, Planck-lOO and Planck-143). 


3.1 Angular power spectrum estimator 

To derive different sets of cross-spectra simulations we use Xpol, a 
pseudo-C^ estimator (PCL) based on the extension to polarization 
of the Xspect algorithm ( [Tristram et al.|2005) . At very low multi¬ 
poles, the PCL estimator is known to be sub-optimal with respect to 
e.g. a Quadratic Maximum Likelihood estimator (QML) (jTegmarkj 
|& de Oliveira-Costa|2001[ |Efstathiou|2006t . This means that the 
variance and correlation of the PCL is expected to be slightly higher 
than for QML resulting in slightly larger distributions for the esti¬ 
mated Cf. However, the implementation of cross-spectra for PCL 
estimators is straightforward and we can easily take into account 
the level of £-hy-{ correlations using Monte Carlo simulations. In 
any case, the definition and validation of the cross-spectra likeli¬ 
hood are independent on the choice of the cross-spectra estimator 
used. 

We produced different sets of Monte-Carlo cross-spectra sim¬ 
ulations. We generated simulated CMB maps on which we add 
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Figure 3. Correlation matrix for the cross-power spectra 100x143 with two different sky coverage: 80% (left) and 50% (right). Each block corresponds to TT, 
EE, BB, TE spectra respectively. { = 0 and ( = \ are not defined and set to zero. 


anisotropic and correlated noise corresponding to four public 
datasets: WMAP (V band), Planck-L¥\ (70 GHz) and two Planck- 
HFI channels (100 GHz and 143 GHz). The CMB signal is con¬ 
structed from the adiabatic ACDM model with cosmological pa¬ 
rameters: Q.i,h^ (the baryon density), Q.ch^ (the dark matter density), 
the amplitude and the spectral index of the primordial power spec¬ 
trum A , and 6 (a parameter proportional to the ratio of the sound 
horizon and the angular diameter distance at recombination), the 
optical depth to reionization parameter t or, equivalently, the red- 
shift of reionization Zre- We also consider primordial tensor modes, 
parametrized by the tensor-to-scalar ratio of the amplitude of the 
primordial spectra r. 

Our reference simulations are generated with a fiducial 
ACDM model based on the Planck 2015 best fit ( [Planck Collab-| 
[oration XIII||2015t with t = 0.078. For the tensor-to-scalar ratio 
we choose r = 0.1 as the fiducial input value. Since it is relevant 
for some validation tests, in particular to check the dependence on 
the fiducial model, we also generated two sets of simulations with 
different input cosmologies: 

model 1: early reionization without tensor modes (Planck 2015 
best-fit with t = 0.09, Zre - 11-2, r - 0) 
model 2: late reionization with high level of tensor (Planck 2015 
best-fit with t = 0.0566, Zre = 8, r = 0.2) 

We estimate the noise angular power spectrum in temperature 
and polarization using the spectra of (I,Q,U) year map differences 
for WMAP and Planck-10. For Planck-HFl, we compute the tem¬ 
perature spectrum from available HFI intensity year map differ¬ 
ences which we rescale according to the number of polarized detec¬ 
tors at each frequency to mimic the polarized noise power spectra. 
The latter ends up very close to what is published in [Planck Collab-| 
[oration VIII[ ( [2015) l. With this procedure, the noise power spectra 
used for the simulations include realistic white noise level and low- 
frequency noise due to systematic and foreground residuals. From 
those power spectra, we derive constrained map realization of noise 
for each simulation. We then scale the noise map by the appropri¬ 
ate relative hit counts in each pixel to simulate the inhomogeneous 
scanning of each dataset. 

We use our pseudo-Cf estimator to produce the six cross¬ 
spectra corresponding to the four datasets: WMAPxPlanck- 
70, WMAPxPlanck-100, WMAPxPlanck-\A3, Planck-lQxPlanck- 


100, Planck-lOxPlanck-143 and Planck-lOOxPlanck-143. For 
each simulation, we construct the TT, EE, BB, TE, TB and EB 
cross-power spectra. The upper and lower panels of Fig. show 
the different noise levels corresponding to the four datasets for the 
E-modes and B-modes spectra, respectively and how the CMB po¬ 
larized power at very low multipoles is directly scaled by the optical 
depth of reionization. 

The plots illustrate the effect of the change form early (zre = 
11.2) to very late (z„ = 8) reionization - which correspond to an 
optical depth of t = 0.09 and 0.566 respectively - for both E and B 
modes below f = 10. In addition, the tensor-to-scalar ratio rescales 
the overall amplitude of the primordial signal in BB. We always 
include the lensing contribution to the B-modes shown as the thin 
line in the lower panel of Fig.|^ 

We do not simulate the impact of foreground contaminations 
in map domain. However, residuals from foreground contamina¬ 
tions are statistically included in our estimate of the noise spectra. 
Moreover, we remove the Galactic plane for the power spectrum 
estimation. We use two sets of Galactic mask based on a threshold 
on the polarized power amplitude of the dust emission and allowing 
for a sky coverage of 80% and 50% respectively. 

The correlation matrices (Fig. are directly derived from 
the Monte Carlo (MC). The level of correlation between multi¬ 
poles depend on the sky cut and the dataset considered. For the 
Planck-lOOxPlanck-143, using 80% sky coverage the correlations 
are weak, lower than 5%, as shown in the left panel of Fig.|^ As we 
will see in the next section (Sect. [3.2] (, for such a large sky coverage 
we can safely neglect the correlations and adapt the full-sky cross¬ 
spectra statistic. For the 50% sky, the correlations are significantly 
higher and can reach the level of 25% (right panel of Fig.[^. 


3.2 Parametrization of the PCL marginals 

The distribution of the PCL estimator is largely non-gaussian in 
the low-f regime we are interested in (see examples in Fig. [Cl[ and 
Fig. [C2[ in the Appendix [C)l and all order moments actually depend 
on the noise powers and on the fiducial model. Leaving aside the 
complicated (and unnecessary) task of defining the full joint p.d.f 
p(C’f), we focus on how to parametrize analytically the individual 
(i.e. marginal) distributions by tweaking the results obtained on the 
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Figure 4. Upper panel: Dependence of the factor on the mask. The 
two distinct groups (full and dashed lines) represents the f^^ values fitted 
to the distributions obtained respectively on the small mask (/sky =0.77) and 
large (/sky =0.49) mask. Lower panel: Dependence of the f^^ factor on the 
fiducial model. The plot shows the f^^ factor as a function of the multipole 
when the early-reionization model (model 1, r = 0.09, solid lines) and the 
late-reionization model (model 2, t = 0.056, dashed lines) are used as input 
cosmology in the simulations. 


full-sky in Sect.j^in the case that a sky cut is applied. We propose 
two different approaches to achieve a satisfactory description. 


3.2.1 Full-sky based approach 

A somewhat heuristic argument used when masking some fraction 
/sky of the sky, is to consider that the ‘number of degrees of free¬ 
dom’ of the associated distribution N = (2£ + 1) is reduced 
asymptotically by the factor ^Hivon et al.||200^ . When the 

mask is apodized by some window, we include the weighting factor 
W 2 /W 4 , where w, is the i-th moment of the weighting scheme, in our 
definition of fs]^y. 

Keeping in mind that cross-spectra do not follow any dis¬ 
tribution and that we are not in the asymptotic regime, we may still 
try to adapt this methodology based on our simulations. We then 
modify our number of modes by N = {21 + keep the gen¬ 

eral full-sky shape of Eq. ij^ and fit for the factor for different 
masks, noise combinations and models. 

The upper panel of Fig.j^shows the factor as a function 
of the multipole t in the case of cross-spectra with different noise 
levels for the small (/s^y =0.77) and larger (/sky =0.49) mask. As 
expected, there is a strong dependency of /(^y® on the mask size. 


The ^ constant in the low-f regime (< 15) and it is 

asymptotically slightly different from the standard /ky factor. This 
can be traced to the fact that we are dealing with polarization that 
involves different Wigner 3j functions than the ones derived from 
temperature. Despite this strong dependence on the mask, the /ky® 
functions derived for the six cross-spectra show a very good con¬ 
sistency for all those different noise levels. 

In the lower panel of Fig.|^we check the model dependency by 
comparing the /Jy® values reconshucted from two simulation sets 
generated with two different input cosmologies: the early reioniza¬ 
tion scenario (model 1, t = 0.09) and the late reionization scenario 
(model 2, T = 0.056). The reconstruction of the /Jy® factor is rea¬ 
sonably stable with respect to the change in the fiducial model used 
and will be considered in the following as independent. 

The stability with respect to the choice of the fiducial model is 
further demonstrated in the Appendix [C|(Fig. [CT| l which shows the 
excellent agreement of the p.d.f s for the Planck-lOOxPlanck-143 
cross-spectrum estimator/or both models while f^^{t) is derived 
from a single one. Note that we choose to display our ‘worse’ case: 
all other cross-spectra parametrisations are even better. 


3.2.2 Edgeworth expansion 

As an alternative to the reconstruction of the /ky®(7’) function, we 
also propose another approach, noticing that only the first three 
central moments contribute essentially to the estimator distribution 
above € > 4. In this case we use the standard (constant) /ky factor 
in A = {21 + I)/sky and, to account for the fact that a fraction of the 
sky is masked, we modify the coefficients of the C*-polynomial 
of the full-sky cumulants (Eqs. |^and|^ to match the ones recon¬ 
structed from the simulations. The new cumulants in this more gen¬ 
eral case take the form: 


K,{cr) = c? 


KdcD = 


1.5(C/)^ + 2C/(Af -H Nf) + NfN^ 

N 


(7) 


6(C?)^ + nCf{Nf + Nf) + lONfNf 

Ki{Ce ) - C( 


Note that this parametrization just depends on constant values 
of the polynomial coefficients. Fig. shows the K 2 (variance) and 
K-i (skewness) reconstructed from the simulations for the two input 
fiducial models (early reionization model in blue and late reion¬ 
ization model in red). The agreement between the full-sky based 
approximation (dashed lines) and the parametrization of Eq. l|^ de¬ 
rived from simulation (solid lines) is excellent. We emphasize that 
the Ki tuning is not mandatory (one may use the one from Eq. ^) 
since it drops rapidly. 

The optimization of the cumulants was performed on a single 
cross-spectrum {Planck-lOOxPlanck-143, model 2). It is however 
robust enough to be used in all other cases as will be demonstrated 
later in the likelihood tests (Sect.j^. 

Now that we have a model-independent parametrization of 
the first cumulants, we proceed in writing an analytical description 
of the estimator p.d.f u sing an Edgeworth Series expansion (e.g. 

Ct-p 


(8) 


Kendall & Stuart 1963 l Using the normalized variable y = 


wEere7r^'^rando^'= s/k^, the truncated expansion reads: 


f{y\Cf,Nf,Nf) = ^(y)(l + ^« 3 (y)), 


where N denotes the normal distribution and is the 3rd order 
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Figure 5. Variance and skewness of the cross-spectra estimators. The plot 
shows the second order cumulants, K 2 (upper red and blue curves) and 
third order k-^, (lower red and blue curves) of the PCL estimator. The blue 
lines correspond to results with the eai'ly reionization scenario (model 1) as 
fiducial model, the red when simulations has the late reionization scenario 
as input cosmology (model 2). The points refer to the cumulants recon¬ 
structed from the Planck-\OidxPlanck-\A?> simulations. The dashed lines 
con'espond to the analytic full-sky derived expressions Eqs. [5|j^ with a 
rescaled N = {2£ + l)/sky factor. The solid lines refer to the parametrization 
of the cumulants based on the Edgeworth expansion defined in Eq. jTJ. 


‘probabilistic’ Hermite polynomial < |Kendall &. Stuart|1963[ t. Each 
cr, K 2 is computed from Eq. ([^ and depends only on Nf,Nf. 

A classical issue when truncating an Edgeworth expansion is 
that, despite being properly normalized to one, it may lead to neg¬ 
ative values so that Eq. is not really a p.d.f and may lead to 
potential problems when constructing with it a log-likelihood func¬ 
tion. We adopt the method proposed by |Rocha et al.|p001| l to al¬ 
leviate this problem. Their idea is based on the solutions of the 
quantum harmonic oscillator, that exhibits non-Gaussianity above 
the ground level. For one extra-level the wave-function (i.e. a p.d.f) 
is of the form: 

P(x) = M(x){oa + ^H 3 (x)^ , (9) 


with tto = yjt - ol - For a mild non-Gaussianity (small a^), which 
is the case in our regime, we expand this equation: 

P(x) = Af(x) Jl -t ^H,ix) + 0(a^)J, (10) 

and equating terms to Eq. we identify: 


Kj 

2 


( 11 ) 


In the following we will refer to ‘Edgeworth expansion’ as this 
regularized form, namely Eq. using Eq. 113 - 

Fig. |C2| shows the agreement between the empirical estima¬ 
tor distributions obtained on simulations and our Edgeworth-based 
parametrization. The agreement is very satisfactory but for the 
C = 2,3 case which would require the use of higher order cumu¬ 
lants. On the other side, introducing some K 4 (kurtosis) term brings 
some oscillations upon all the multipoles. This would not be desir¬ 
able since the first two accessible multipoles have generally a very 
low SNR due to 1// noise and large cosmic-variance and can be 
disregarded without a sizable loss of information. 


4 CROSS SPECTRA-BASED LIKELIHOODS 

With these tools in hand we now proceed in constructing the likeli¬ 
hood of a given model, which means inverting the (unknown) joint 
and possibly multi-field PCL estimator distribution given the true 
value Cf. We first discuss the simple but frequent single-field case 
with a small mask for which we give analytical formulas. We then 
define a more general solution based on the modification of the 
H&L approximation to construct a general likelihood solution for 
the combination of the temperature and polarization fields account¬ 
ing for correlations. 


4.1 Single field approximations neglecting correlations 

As a first solution, we can build our real-case likelihood from the 
parametrization of the marginalized estimator distribution p{Ce) 
described in Sect. |3.2| This approximation is accurate when the 
masked sky-fraction is low (typically below 20%) so that the £- 
hy-£ correlations can be safely neglected. The likelihood function 
is defined as the product of the probability density functions 
(cfr. Eq. 1^): 


t-max 

y^^\Cf\Ct,N^,Nf)= f] pf\Ct), (12) 

t=tmm 

where the Q represent the values measured on data. The p^^ func¬ 
tions are implicitly dependent on N^, Nf, Cf and they can be 
defined according to the chosen analytical parametrization as de¬ 
scribed in Sect. l3.2.Tl and Sect. 13.2.^ 

Note that this approximation derived from the full sky formal¬ 
ism is easily defined for a single field, that is to say when the cross¬ 
spectra are derived from the combination of the same temperature 
or polarization field, e.g. the E-modes cross-spectra. A combined 
analysis that accounts for all the temperature and polarization fields 
is difficult to define analytically since correlations between different 
fields (TE, TB, EB) cannot be neglected and higher order moments 
of the Cl distribution must be accounted, making the analytical so¬ 
lution difficult to handle in this more general case. 


4.2 General multi-field approximation 


Here we present a more general formalism to define a cross-spectra 
likelihood for the analysis of the CMB data at large angular scales 
that allows to deal with realistic cases of incomplete sky coverage 
taking into account the t - £ correlations. This likelihood can also 
be easily generalized to a multi-fields likelihood that combines the 
temperature and polarization fields T, E and B, accounting for the 
field-field correlations TE, TB and EB. 

In order to model the non-Gaussianity of the Q estimators, 
the approximation that we propose is based on the modification of 
the Hamimeche&Lewis likelihood (H&L) ( [Hamimeche & Lewis| 
2008 1 , adapted to work for the cross-spectra and at low-f. 


The general form of the H&L likelihood is defined for auto¬ 
spectra at intermediate and small scales {£ > 30) ( [Hamimeche &| 
|Lewis|2008| (: 

- 2ln,^(Cf\Ci) = (13) 


The is the inverse of the Q-covariance matrix that allows 

to quantify the £ - £ and the correlations of the T, E, B fields. The 
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vector [Xg\i is the H&L transformed Q vector defined as: 


[X,]t = vecp(c}/^U(g[D(P)])UTc};,"). 


Ineq0C}(,^ 


is the square root of the Q matrix: 


(14) 



' f-^TT 

r-TE 


Q = 

f-TE 

f-'EE 

f-'EB 


f'TB 

\yi 

f'EB 

r'BB 

) 


(15) 


for a given fiducial model and the function (g[D(P)]) refers to the 
transformation: 

g(x) = sign(x - 1) -^(2(v - ln(x) - 1)), (16) 

applied to the eigenvalues of the matrix P = where 

Cmad and Cdata are, respectively, the matrices of the sampled Cf 
and the data. This approximation has been shown to be robust with 
respect to the choice of the fiducial model ( [Hamimeche & Lewis| 
|2008| >. The problem is that for cross-spectra and at large angular 
scales the P matrix is no longer guaranteed to be positive definite. 
In fact, as shown in Eq. Q that describes the distribution of the 
cross-spectra estimators Cc, the Cdata can be negative. In order to 
solve for this issue, we propose a modification of the H&L like¬ 
lihood that consists in adding an effective offset 0 ( to the cross¬ 
spectra. This mimics the noise bias of the auto-spectra and makes 
the offset-cross-spectra distribution very similar to the auto-spectra 
distribution used in the H&L approximation. We re-define each Cf 
matrix (eq. [Tg as: 



'ey + of 

f-TE 

f-'TB 

'-'£ 


cf ® ^ 0(Cf ®) = 

r'TE 

cy + oy 

f-'EB 

'-'£ 

(17) 


/^TB 

\ 

f-'EB 

^£ 

qBB , BB 



so that: 


[X,(Cf'^)]f ^ [OX,]e = [X,(0(Cf'^))]f. (18) 

The new offset H&L likelihood (oHL hereafter) reads: 

- 2ln^{Ci\C^^'^) = Y}OX,]][M;f^]u'[OX,]e'. (19) 

w 

The variable transformation g(x) is now modified for the cross¬ 
spectra to regularize the likelihood around zero so that Eq. GD 
now reads: 


g(x) -> sign(x)g(lxl). (20) 

The offset function of^ can be derived from simulations. We 
estimate the offsets from the MC distributions ensuring that the P 
matrix reconstructed is positive definite for more than 99% of our 
simulations. Given that the offsets are needed to shift the Cf distri¬ 
butions for each field T, E, B to avoid negative eigenvalues on the 
P matrix, the offset functions depend on the shape of the Cf distri¬ 
bution at each {. In particular, the offsets depend on the noise levels 
of the maps involved in the cross-spectra and on the mask used. 
In fact, the tails of the Cf distributions at each t are more negative 
when the noise is higher and when a larger mask is applied. The 
plot in Fig. 1^ shows how the offset functions change for different 
combinations of noise levels in the case of the six cross-spectra 
considered: from the highest of the WMAPxPlanck-70 in orange 
to the Planck-lOOxPlanck-143 in green. 

Moreover, the offsets also depend on the fiducial model, as. 



Figure 6. The offset functions of the oHL likelihood for the E-modes 
cross-spectra and for the six different combinations of noise levels. 



Figure 7. The offset functions of the oHL likelihood for the Planck- 
HFlxPiafick-143 B-modes cross-spectra (solid) and the E-modes cross¬ 
spectra (dashed). The different colors refers to the different fiducial mod¬ 
els used in the simulations: black is the early reionization scenario without 
tensor modes (model 1), red is the late reionization scenario with tensors 
(model 2) and blue is the Planck 2015 best fit with r = 0.1. 


in general, an higher signal-to-noise implies that the C; distribu¬ 
tions have a smaller shift to negative values. Fig.[7]shows the offset 
functions obtained from simulations generated with different fidu¬ 
cial models for the E-modes (dashed) and B-modes (solid) for the 
Planck-lOOxPlanck-143 cross-spectra. The black lines refers to the 
early reionization scenario without tensor modes (model 1), the red 
lines to the late reionization scenario with tensors (model 2) and 
the blue lines to the Planck 2015 best fit with r = 0.1. The shape 
of the offsets is consistent for the three different cases and, given 
the very different levels of signal considered, the dependence on 
the fiducial model is mild. As we will show in Sect. |5] and Sect. 
the likelihood results on the cosmological parameters reconstruc¬ 
tion are robust with respect to the choice of the fiducial model used 
to define the offset functions. 


The in Eq. 1 19i is the inverse of the cross-spectra 


C^^"-covariance matrix that can be estimated for a given theoret¬ 
ical fiducial model through Monte Carlo simulations such 

that: 




= <((Cf),, 


^XY fid\/ 

■ 'aim 


^XYfid\. 

j/MC-> 


( 21 ) 


where Cp = and X,Y = {T,E,B]. 
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Since it will be useful in the following, we also report the 
equations of the modified oHL likelihood in the case of the single 
field approximation. In particular, we are interested in applying the 
method to the polarization EE-only cross-spectra 
for which the oHL likelihood is defined by 

- 21n,if = ^ [OX,]f [Mf (22) 

ef' 

where: 

[x.r ^ [oz,]f = 

(23) 

and: 

0(Cf) = (Cf -H Of). (24) 

C^Efid ^ ^ee QEEmod respectively, the spectra of the fiducial 
model, the data and the variable spectra for the likelihood sam¬ 
pling, and is the effective offset. Also, the covariance matrix to 
account for the multipole coupling in this case is defined by: 

[Mf]u' = - C,“^'‘')((C“)^-f - Cf (25) 

5 SINGLE FIELD RESULTS 

We first present the results in the case of the single field approx¬ 
imation. As single field we choose the E polarization and we 
build the EE cross-spectra likelihoods to constrain the t parame¬ 
ter, since it is relevant for the analysis of present and future CMB 
data. We construct the three different single field cross-spectra like¬ 
lihoods derived from the formulas in Sect. the general ana¬ 
lytical parametrization derived from full-sky based approach, the 
parametrization based on the Edgeworth expansion approximation 
to describe the cumulants of the cross-spectra distribution and the 
oHL single-field likelihood. 

In order to compare the three methods, we focus on the small 
sky cut case, where the cross-spectra simulations are generated by 
applying a mask with /sky =0.8. The t-by-t correlations are weak 
and the analytic approximations are reliable. This comparison is 
useful not only as a validation test of the different methods but 
also to demonstrate that correlations can indeed be neglected in 
the parametric case. To construct the single field oHL cross-spectra 
likelihood we use Eq. \22) where the £-C correlations are encoded in 
the cross-spectra covariance matrix of Eq. ( |25^ . Eor each of the six 
cross-spectra considered, the covariance matrix is computed from 
the Monte Carlo average of 10.000 E-modes simulations generated 
with a fiducial input cosmology corresponding to the Planck best-fit 
2015 with r = 0.078 and tensor modes with r = 0.1. We estimate 
the offsets from our reference simulations as described in 

Sect.lL^andEig.l^ 

Note that, as pointed out in Sect. |3.2| and Sect. |4.2| the 
parametrization used to define the analytical approximations and 
the definition of the offset functions of the oHL likelihood are both 
robust with respect to significative changes of the t parameter in 
the fiducial model (Ar^y = 0.03 » fr^). However, in general, 
as it is the case for the HL likelihood analysis at smaller scales 
( [Hamimeche & Lewis|2008) ), it is a good choice to use a fiducial 
model close to the "true" model to compute the covariance matrix. 

The likelihood sampling is done by computing the 
with the CAMB cod^ varying t in the range [0.01,0.15] with 

' http://camb.info/ 


a step At = 0.001, fixing the other parameters to their Planck 
2015 best-fit values and rescaling . The degeneracy between 
T and the scalar amplitude parameter A, is in fact broken by fix¬ 
ing accordingly the amplitude of the first peak of the TT spectrum 
Att = at f = 200. More general results based on joint con¬ 

straints of the T and As are presented in Sect. |6.2.7| 

To compare the three likelihoods, we choose events (i.e. one 
C( vector sample) at random from the set of Planck-lOOxPlanck- 
143 simulations and construct for each t independently the 
marginal likelihoods with the three different methods, setting each 
time all Cf values other than this multipole to their true values. Eig. 
[^displays a typical case. Here are some comments that we derive 
from the observation of many samples: even-though the sample Ct 
may get negative values, due to noise and low signal, the likeli¬ 
hood of any negative true power value is unphysical and is equal 
to 0. This case does not happen in practice since in cosmological 
parameter estimation the Boltzmann code always propose positive 
spectra. The Edgeworth-based method shows some oscillation for 
the very first multipoles, generally for £ = 2,3. This is due to the 
very steep raising of the distributions at the very beginning (see Eig 
|Cl[ l which leads to some small negative ‘ringing’ effect in the trun¬ 
cated expansion. The method introduced in Sect. |3.2.2| mitigates 
the effect but does not completely cure for it and regarding this 
aspect the full-sky based method and the oHL likelihood gives a 
better approximation. Overall Fig.j^shows the excellent agreement 
among the three likelihood methods in recovering the minimum 

= —2\n[^{C'^)IJ^max\ for each multipoles with comparable 
accuracy. 

As a further validation test, we check the bias of the likelihood 
against our set of 10.000 Monte-Carlo simulations. For each sim¬ 
ulation, we derive the distribution of the maximum likelihood of t 
for £ < 20. For the full sky based likelihood and for the Edgeworth 
expansion likelihood we remove multipoles 2 and 3 -which do not 
carry much information due to the cosmic variance level- since their 
p.d.f parametrization is less accurate, as shown in Sect. |3.2| For 
oHL we consider £ = [2,20]. 

Figure shows the distribution of the maximum probability 
over the Monte Carlo simulations for the full sky based likelihood, 
the Edgeworth expansion likelihood and the oHL likelihood. All 
three approximations recover the input value Tyy = 0.078 used to 
generate our reference simulations, showing that the three likeli¬ 
hoods are unbiased. In Table [T] are reported the best fit values and 
the error bars on the estimation of the t parameter. The error bars 
are computed as the standard deviation of the maximum probability 
f for the three likelihoods. Since the oHL likelihood accounts for 
the £-by-£ correlations while the full-sky based likelihood and the 
Edgeworth expansion approximation do not, the fact that the three 
methods give compatible results in terms of error bars confirm that 
the level of multipole correlations for a small sky cut is low and 
does not have an impact in the reconstruction of the t parameter. 
Note however that both the analytical approximations are slightly 
sub-optimal with respect to the oHL likelihood, by a factor of = 4% 
for the full-sky based likelihood and by a factor of = 7% for the 
Edgeworth expansion likelihood. These results hold in general for 
all the cross-spectra considered. 

Finally, it is useful to assess the stability of the results obtained 
with the oHL likelihood with respect to choice of the offset term. 
Indeed, changing the offsets both could bias the peak of the pos¬ 
terior distribution and change its width. As described in Eq. j23| l, 
the offset ensure the H&L transformation to be definite and too 
small otfsets may leak to undefined likelihood. On the opposite, 
a overestimation of the otfset value has limited effect on the peak 
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Figure 8. Comparisons of the three likelihoods methods developed for the low-correlations case (/sky = 0.8): red is for the Edgeworth expansion method, blue 
for the full-sky based method and green for the modified Hamimeche-Lewis approximation (oHL). The arrow represents where our random sample has fallen. 
Due to noise, it can be negative. The ordinate is Ax^ = -21n[=5f(C*)/oSfmox]- 


Table 1. Comparison of the best fit estimation of the reionization optical depth f and error bars CTf for the three likelihood methods from simulations (see also 
Fig.[9|. The errors are computed as the standard deviation of the maximum probability f over a set of 2000 simulations. The input value used in the simulations 
is = 0.078 and /^ky = 0.8. 


Cross-spectra 

(f ± o-f 

(f ± o-f)''''*'’"'®"* 

(f ± 

WMAPxPlanck-70 

0.0777 ±0.0116 

0.0768 ±0.0121 

0.0774 ±0.0110 

WMAPxPlanck-m 

0.0777 ± 0.0088 

0.0768 ± 0.0092 

0.0773 ± 0.0086 

WMAPxPlanck-m 

.00774 ± 0.0086 

0.0765 ± 0.0089 

0.0772 ± 0.0084 

Planck-1 dxPlanck- 100 

0.0776 ± 0.0077 

0.0773 ± 0.0079 

0.0774 ± 0.0074 

Planck-1 dxPlanck- 143 

0.0775 ± 0.0075 

0.0771 ±0.0078 

0.0774 ±0.0071 

Planck- 1 OOxPlanck-143 

0.0777 ± 0.0054 

0.0778 ± 0.0055 

0.0781 ±0.0051 


distribution. Figure^] shows that the impact of a factor of two in 
the estimation of the offsets amplitude is negligible on the posterior 
distribution of the r parameter. The figure illustrates two represen¬ 
tative cases of r posteriors obtained with the highest and lowest 
noise configuration from our simulations. Note that a change in the 
offset of this type could arise if the fiducial model used to derive 
the offsets is very different from the best fit model, as illustrated in 
Fig. § The fact that this change has practically no effects on the 
posterior distributions demonstrates that the definition of the oHL 
likelihood is robust with respect to the offset reconstruction. Also, 
the effect on the width of a change in the offset is very weak, mean¬ 
ing that the offsets, as expected, do not affect the estimation of the 


error bars. The same results hold true for all the cross-spectra con¬ 
sidered. In general, our offset terms are well defined and the oHL 
results are robust with respect to the offset choice. 


For a smaller sky coverage, the three likelihoods remain un¬ 
biased even if the C-by-£ correlations get larger. In this case, the 
parametric likelihoods are less optimal: their variance increase by 
about 10% compared to the oHL likelihood for the 50% mask. 
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Figure 10. The plots show the effect of changing the offsets amplitude by a 
factor of two on the posterior distribution of the r parameter. The top panel 
refers to the WMAPxPlanck-10 cross-spectra, while the bottom panel to the 
Planck -1 OOxPlanck -143 cross-spectra, 


cross-spectrum we construct the full [T,E,B] covariance matrix of 
Eq. by computing the Monte Carlo average of 10.000 sim¬ 
ulations generated with a fiducial input cosmology corresponding 
to our baseline Planck 2015 best fit with t = 0.078 and r = 0.1. 
The offsets functions are derived from the same simulations as de¬ 
scribed in Sect. |4.2| For each cross-spectra we therefore add the 
offsets oY, Of Of^ to the diagonal elements of the Ce matrix as 
defined in Eq. tni- 


6.1 Constraints on r 


Figure 9. Distribution of the maximum probability for the analytic full- 
sky based likelihood {top), for the analytical parametrization based on the 
Edgeworth expansion (middle) and for the oHL likelihood (bottom) on each 
E-modes cross-spectra on 80% of the sky 


6 RESULTS FOR CORRELATED FIELDS 

This section is dedicated to the results obtained with the full tem¬ 
perature and polarization oHL likelihood (Eq. ([T^). One of the 
main advantages of the oHL method relies in fact on the possi¬ 
bility to include in the analysis both the correlations between the 
[T,E,B] fields and the multipole correlations. Since the simulations 
used for each cross-spectrum are built with realistic noise levels 
as described in Sect.[^ the forecasted estimates on the t, r and Aj 
parameters from the low-^ analysis presented here are realistic for 
current CMB experiments. 

We build the {TT, EE, BB, TE, TB, EB) oHL like¬ 
lihood for the six different cross-spectra: WMAPxPlanck- 
70, WMAPxPlanck-100, WMAPxPlanck-143, Planck-lOxPlanck- 
100, Planck-1 OxPlanck-1 A3, Planck-100xPlanck-lA3. For each 


Firstly, we study the impact of including the T,E,B cross-spectra 
and their correlations on the estimation of the optical depth to reion¬ 
ization T, compared to the single-field EE analysis described in the 
previous section Sect.|^ The sky fraction is fsky = 0.8 and the mul¬ 
tipole range used is t = [2,20]. As shown in Fig. [TTjand Table 1^ 
the combined analysis gives unbiased results on the estimation of r. 
As expected, adding the temperature and the tensor modes and all 
the possible correlations gives results very close to the single-field 
EE analysis since the relevant physical information related to t is 
essentially encoded in the EE-spectra. However, the full tempera¬ 
ture and polarization analysis leads to a slight improvement in the 
estimation of the t error bars, which is of about a few percent for 
the Planck-l00xPlanck-lA3 analysis. 

We compare the t posterior distribution of the oHL likelihood 
to the one from the pixel-based likelihood. We implement a pixel- 
based likelihood by using a combination of the maps at Planck-100 
and Planck-l A3 from the same simulation set that we used to gener¬ 
ate the Planck-HFl cross-spectra. Both methods are therefore based 
on simulations with the same noise characterization. A typical case 
is given in Fig.[T^ As expected, the oHL likelihood approximation 
is slightly sub-optimal with respect to the pixel-based likelihood 
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Figure 11. Validation of the oHL multi-fields likelihood. The plots show 
that the oHL likelihood computed combining the T, E and B fields and ac¬ 
counting for both multipole and fields correlafions gives unbiased results 
on the estimation of the optical depth to reionization parameter t. The top 
panel shows the t posterior for the six different cross-spectra when 20% of 
the sky is masked (fsky = 0.8), while the bottom panel shows the results 
for a bigger mask with fsky = 0.5. The dashed line refers to the input value 
tfid = 0.078 used in the simulations. 


Table 2. Results on the estimation of the r parameter with the full tem¬ 
perature and polarization oHL likelihood. The fiducial model used in the 
simulation is the Planck 2015 ACDM best fit with Tfjj = 0.078. The table 
shows the comparison between the r estimates (best fit f and en'or bars CTf) 
obtained with two set of simulations with different sky cuts: the small mask 
with with fsky = 0.8 and a bigger mask with fsky = 0.5. 


Cross-spectra 

f ± O-f {fsky = 0.8) 

f ± o-f (fsky = 0.5) 

WMAPx/’/anck-70 

0.0750 ± 0.0108 

0.0761 ± 0.0203 

WMAPxPlanek-mO 

0.0769 ± 0.0075 

0.0764 ± 0.0121 

WMAPxP/anck-143 

0.0769 ± 0.0079 

0.0770 ±0.0116 

Planck-lOxPlanck- 100 

0.0783 ± 0.0069 

0.0776 ± 0.0105 

Planck-lOxPlanck- 143 

0.0784 ± 0.0065 

0.0763 ± 0.0101 

Planck- 1 OOxPlanck- 143 

0.0780 ± 0.0049 

0.0788 ± 0.0069 


which is not an approximation and is build to be statistically opti¬ 
mal. Note however that the error bars obtained with the oHL like¬ 
lihood are comparable with the optimal estimate obtained by using 
the pixel-based approach at better than 15%. 

Finally, we use the combined oFIL likelihood to test the re¬ 
sults with a different sky cut. We consider a severe cut at 50% 
(/itv = 0.5). This is a more complicate case to deal with since 
the £-hy-{ correlations are stronger. Also, the shape of the distri¬ 
butions of the C{ estimators at each { is affected by the smaller 
sky coverage, leading to more negative tails. We generate the offset 
functions for each cross-spectra as described in Sect. |4.2[ using our 
reference simulations masked at 50%. The results are summarized 
in tablej^and in the bottom panel of Fig.[TT]that shows the t pos¬ 
teriors for each of the six cross-spectra. Even in this more complex 
case, the oHL likelihood analysis is unbiased. As expected, since 
we are considering a smaller sky fraction and non-negligible mul¬ 
tipole correlations, we recover bigger error bars with respect to the 
fsky = 0.8 analysis, with a degradation of =: 30% for the Planck- 
lOOxWauck-143. 



Figure 12. Comparison of the posterior distributions of the r parameter ob¬ 
tained with the full temperature and polarization oHL likelihood (blue) and 
with the pixel-based likelihood (green). The plot shows a typical example 
from the Planck-HVl simulation set. 


6.2 Joint estimation of t, r and Aj 

Using the full combined analysis, we can construct multi¬ 
dimensional constraints on parameters. In particular, we focus on 
the correlations between the optical depth and the amplitude of 
the scalar fluctuations A, and between the optical depth and the 
tensor-to-scalar ratio r which are relevant for the future analysis of 
CMB data at large angular scales from e.g. Planck. In both cases we 
perform the full analysis using the Planck-100xPlanck-lfi3 spec¬ 
tra which corresponds to the lowest noise frequency combination 
and it can be used to make realistic forecasts for current and future 
CMB experiments. We consider a sky cut with fsky = 0.8. 

6.2.1 Joint estimation of T and As 

Using the temperature power spectrum only. A., and t are strongly 
degenerated. Indeed, the amplitude of the first acoustic peak of the 
CMB temperature power spectrum directly measures Att = 

Using polarization data at large angular scale helps breaking this 
degeneracy. So far we fixed the degeneracy direction by rescaling 
the temperature spectrum, fixing A^j, accordingly to the variation 
of T in the likelihood sampling. Here we let A, free to vary. The 
results from the simulations, using the Planck-lOOxPlanck- 143 full 
oHL likelihood are summarized in the left panel of Fig. [T^ The 
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Figure 13. 2D-distribution of the maximum likelihood for T-r (right panel) and r-A., (left panel). The plots show the joint constraints obtained with the full 
temperature and polarization oHL likelihood on 2000 simulations of the Planck-lOOxPlanck-143 cross-spectra. The fiducial input parameters used in the 
simulations are: Ty/j = 0.078, ry/j = 0.1 and [/n(10*®Aj)]y,rf = 3.09. 


plot shows the 2D histogram of the best fit values for the whole 
set of simulations in the t-A, projection. The full oHL likelihood 
correctly recovers the inputs values for t and as well as error 
bars compatible with the MC dispersion. 


6.2.2 Joint estimation ofr and r 

The CMB power spectra at large angular scales, in particular the 
E and B polarization modes, are affected by how the reioniza¬ 
tion process proceeded and lasted. Thus, as shown in Fig. the 
power at large scales (low-f) in the B-modes spectrum is directly 
related to the reionization optical depth. Indeed, the amplitude of 
the B-modes spectrum reionization bump scales with t“: 

T-Cf^ 2 o- '^he amplitude of the B-modes spectrum of course also 
depends on the amount of the primordial tensor perturbations, there 
is a degeneracy between the t and r. It is therefore interesting to de¬ 
rive joint estimates of these parameters. 

We compute the joint T-r constraints with the full oHL likeli¬ 
hood on the set of 2000 simulations of the Planck-lOOxPlanck-143 
cross-spectra with an input cosmology corresponding to the Planck 
2015 best fit for the base ACDM parameters with t = 0.078 and 
a tensor-to-scalar ratio of r = 0.1. The multipole range used is, as 
usual, { = [2,20]. The results of the oHL likelihood sampling on 
simulations are summarized in the right panel of Fig. The plot 
shows the posterior in the T-r plane from the oHL and from which 
we can see that the oHL likelihood correctly recovers the param¬ 
eters T and r compatible with the input values used in the simu¬ 
lations. As regarding the error bars, the forecasted Icr error for t 
in the case of the highest resolution channels of a Planck-\fk& ex¬ 
periment is a-™xi 43 _ Q Q 05 J Pqj- tensor-to-scalar ratio in the 
multipole range considered, we find a-Jooxi 43 _ q Q 9 Note that in 
our analysis, we consider a correlated noise model. This noise char¬ 
acterization, which is more realistic with respect to a simpler white 
noise modeling, implies a rising of the noise level at low multipoles 
due to the 1// noise correlations (see Fig.[^. Therefore, in particu¬ 
lar in the case of a low signal scenario, the correlated noise at large 
scales can eventually dominate over the cosmic variance inducing 
a worsening of the constraining power proportional on how steep is 
the rising of the correlated noise at low multipoles. 


7 CONCLUSIONS 

In this paper we presented a new approach for the analysis of 
the CMB polarization data at large angular scales based on cross¬ 
correlation in spectra domain. Using cross-spectra with respect to 
the auto-spectra and, in general, to the pixel based approach used 
so far in the CMB analysis at large angular scales has many advan¬ 
tages, in particular in the case of a realistic CMB experiment that 
accounts for anisotropic noise and a sky cut needed to minimize the 
foreground contamination. In fact, by using cross-frequency/cross- 
dataset CMB spectra, the noise biases and the systematics specific 
to a given frequency/dataset are removed. Also, the possible fore¬ 
ground residuals can be minimized and the information encoded in 
different frequencies/datasets can be combined efficiently. 

The cross-spectra estimators are non-Gaussian at low multi¬ 
poles especially when applied on cut-sky. We generalized the ap¬ 
proximation made in |Hamimeche & Lewis|p008t to accommodate 
for cross-spectra at very low multipoles. This likelihood (oHL) can 
easily handle the correlation between CMB modes (TT, EE, BB, 
TE as well as TB and EB) and between multipoles and gives er¬ 
ror bars less than 15% larger than the optimal pixel-based method. 
The oHL likelihood shares the same robustness with respect to the 
choice of the fiducial model as the H&L approximation (see dis¬ 
cussion in |Hamimeche & Lewis] {2008[ l). We compared the oHL 
likelihood to the analytical parametrization of the estimator distri¬ 
bution which can be used as a quick likelihood solution in the case 
of a single field analysis with small sky cuts so that correlations can 
be safely neglected. 

We generated different sets of simulations that we used to 
construct and validate the likelihoods, proving that all the meth¬ 
ods are unbiased and can accurately constrain the optical depth 
to reionization parameter t. Also, we showed that the oHL like¬ 
lihood gives accurate constraints of the joint estimation of the t 
parameter, the tensor-to-scalar ratio parameter r and the amplitude 
of the primordial scalar perturbations A„. Our simulations account 
for anisotropic correlated noise, beam, mask with the characteris¬ 
tic of a realistic CMB experiment as WMAP and Planck. In order to 
validate our likelihoods for different noise levels, we generated sim¬ 
ulations for cross-frequency spectra with different resolution, from 
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the lowest, WMAPxPlanck-10, to highest, i.e. Planck-lOOxPlanck- 
143. 

Optimal foreground cleaning is beyond the scope of this paper 
but foreground residuals, in particular synchrotron and dust, must 
be quantified in a realistic CMB analysis. In this paper we work 
with cleaned CMB maps but we account and propagate the uncer¬ 
tainties related to the foregrounds removal by using in our simu¬ 
lations realistic estimates derived from public data. The correlated 
noise term that we include in the simulations in fact is drown from 
real data and can be taken as a good proxy for a realistic combina¬ 
tion of noise, systematics and foregrounds residuals, in particular 
at low multipoles. 

The cross-spectra likelihood approach presented in this paper 
is a powerful and efficient tool for the analysis of the CMB data at 
large angular scales. It allows to minimize the impact of the experi¬ 
mental residual systematics (from both instruments and foreground 
contamination) while providing nearly-optimal constraints on the 
estimation of the t, r and cosmological parameters. 
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APPENDIX A: CROSS-SPECTRA DISTRIBUTION ON 
THE EULL SKY 

We consider the product of a pair of correlated central gaussian 
random variables: 


x = ZaX Zb (Za, Zb) ~ 0, V) 


(Al) 


where n is a generic vector and the covariance matrix is written in 
the standard form: 


V = 


cr; 

pcA<r B 


^ per ac B 


(A2) 


Standard probability rules allow to compute its p.d.f jGr- 
|ishchuk|1996^ : 

px 


fAx) = 


1 


_g(l -p^)o-ACrBXa 


\x\ 


no- actB Vl -P^ 
whose characteristic function (Eourier transform is): 

<f>At) = E\e-"]= ^ 


(1 -fA)o-AcrB 


yjl - lipa-AO-Bt + (1 -p^)(T\crlfi 


(A3) 


(A4) 


The sum of N such independent variables X = 2™ j x, has 
therefore the characteristic function: 


0 x(f) 


[ 1 - lipo-AO-Bt -I- (1 - p^)cri(Tl 


(1 - p^)(Tlcrl(t- 


^)] 


-N/2 


-)(t + 


(l-p)crACrB {l+p)crACrB 
To obtain the X p.d.f we inverse-Rourier it: 


fx(x) 

In 


I 




(A5) 

-Njl 


(A6) 


■I 


(f- 


-)(t + 


(1 -p)CAa-B 
and perform the change of variable t ^ 1 + 
px 


(1 +p)crAcrB 

ip 


Nil 


dt. 


(1 -p^)<rA(TB 


to obtain 


fx{x) cc P^)caCb J' 


P + 


1 


(1 -p^)crAcrl 


-Nil 


dt. (A7) 


Then making use of the Basset integral ( |01ver et al.| |2010| 
Eq. 10.32.11) and reintroducing the normalization, we get: 


px 


|^|(iV-l)/2g(l-p2)o-^0-B^^ 


(<V-l)/2 


fx(x) =- 


\x\ 


(1 -p^)a-A<TB 


2 (iv-i )/2 yp:T(NI2) xjl - p2(o-aO-b)<"+‘W, 


(A8) 


where R refers to the gamma function and is the modified Bessel 
function of second kind and order v = (N - l)/2. We can check a- 
posteriori that we recover indeed Eq. jA^ for the N = I case, 
which justifies Eq. l |A4[ l. 

We now have all in hands to consider a full-sky Ax B cross- 
1 ^ 

spectrum Q = ——- ^ where for (isotropic) noise 

l=-m 

power the covariance matrix reads 


V = 


Cf + Nf 

^th 


/~'\h 


(A9) 
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Its p.d.f for a given I therefore reads: 


MCe) = NfxiNCt) (AlO) 


where N = 2t + \, and, in Eq. I |A8^ : 


<Ts=^Cf + N^ 

p- , ' 

^(Cf + NfXCf + Nf) 


(All) 


This formula is similar to the one given (but not derived) in |Perciv^ 
|& Brown] ( |2006| Eq.l9) for the TE distribution. 

The characteristic function of the cross-spectrum estimator is: 


0 c(f) = 


(A12) 


APPENDIX B: AUTO-SPECTRA AND CROSS-SPECTRA 
STATISTICS COMPARISON 



Figure Bl. Ratio of the variance of the A - B coss spectrum estimator to 
the one from the auto-spectrum of the optimally combined map. The colors 
indicate different noise combinations according to the following scheme: 
WMAPxPlanck-10= orange, WMAPxPlanck-lOO = gold, WMAPxPlanck- 
143 = purple, Planck-1 OxPlanck-WQ = red, Planck-lOxPlanck-143 = blue, 
Planck-lOOxPlanck-143 = green. The dashed line recalls that most of the 
interesting information in EE about reinization is contained below f < 15. 


It is instructive to study the respective merits of the auto and cross 
spectra approaches to estimate a single field power-spectrum. We 
concentrate here on the full sky case where the auto and cross spec¬ 
tra estimators can be handled analytically. The main conclusions 
hold essentially on a cut-sky too. 

Let us first recall some properties of power-spectrum estima¬ 
tion using auto-spectra. We consider the measurement of a Gaus¬ 
sian field (of power-spectrum Cf ) over the full sky with an instru¬ 
ment which has an isotropic noise (of power-spectrum N'f) uncor¬ 
related to the signal. 

According to the spectral theorem, the decomposition of the 
map onto the (orthogonal) spherical harmonics basis yields a set 
of independent Gaussian random variables: a/,„’s. Eor a given mul¬ 
tipole (, this variable x = fl/„, follows a Gaussian distribution of 
null mean and variance cr^ = C* + Ne. At a given { , from a set 
of A = 2^ -I- 1 measured harmonic coefficients {xj the Maximum 
Likelihood estimator of C* is the empirical variance: 

(Bl) 

/=1 

This estimator’s distribution can be computed analytically by 
noticing that since x, follows Aff.; 0, C* -l- Nf) it can be written as 

, Cf + Af 

(B2) 

where y follows now a distribution. Then from the analytic x- 
p.d.f and standard probability transformation rules one can obtain 
analytically the full Q p.d.f (which is a T one).|^ 

We do not need any elaborate expression to compute the first 
order moments. Instead we use the scaling property of the cumu- 
lants (which are equivalent to central moments up to the third order) 
which states that if X has some cumulants Kj(X], the cumulants of 
a linear transformation Y = aX + b are a:,(T) = a‘Ki(X) -t- bd-. S ince 
the first cumulants of &x\i distribution are ^ = A(L 2,8...), Eq. 02| 


^ We emphasize we are dealing for the moment with the estimator distri¬ 
bution, not the posterior or likelihood one. 


gives immediately: 

Ki(Ce) = E [Ce] = A - Nf = C? 


K 2 {Ct) = Var[Q] = 


A 

C?+Af 


A 


2N-- 


2(Cf + Nff 


A 


(B3) 

(B4) 


K2{Ce) = E[{Ci-Cy] 


A 


8 A = 


8 (C'‘’ + N'tf 




(B5) 


Any error on the noise level estimation A^ bias accordingly 
the estimator (/tj) and we can recognize the cosmic variance ex¬ 
pression {K 2 ). 

The clear advantage of using cross-spectra has already been 
emphasized: the estimator is unbiased whatever knowledge we 
have of the noise spectra, cf. Eq. l[§. 

However in order to investigate its discriminating power, we 
also need to consider its variance, and to a lesser extent its higher 
order moments. What do we loose statistically using cross-spectra 
over using an auto-spectrum assuming a perfect knowledge of the 
noise? The answer depends on the relative levels of the signal (C*) 
and noise spectra ) and on the t range under considera¬ 

tion. To get some further insight, we consider the EE field from the 
Planck 2015 best-fit ^Planck Collaboration XIII|2015[ ) and variety 
of realistic maps noise-levels. We focused on three particular pub¬ 
lic datasets: WMAP (V band), Planck-10 (70GHz) and two Planck- 
HEI channels (lOOGHz and 143GHz). 

We then consider the variance of the cross-spectrum estimator 
Eq. (0 and compare it to the one obtained on the auto-spectrum of 
an optimally inverse-variance combined map, i.e. with noise 


1 


1 1 


(B6) 

"f ”f 

We show the ratio of these quantities for the different noise- 
pair combinations on Fig. |Bl| 

The variance increase is important in the case of two very dif¬ 
ferent noise levels (for instance WMAPxPlanck-Y{F\) and moderate 
when they are similar (WMAPxPlanck-10, Planck-lOOxPlanck- 
143). This can be understood from the variance formulas Eqs. l |B4l > 


© 0000 RAS, MNRAS 000, 000-000 




























Large-scale CMB cross-spectra likelihoods 15 

APPENDIX C: p.d.f PARAMETRIZATION 

In this section we show the excellent agreement that is obtained 
when comparing the parametrization of the EE PCL estimator dis¬ 
tribution defined with the full sky based approach and the Edge- 
worth expansion method described in Sect. |3.2.1| and Sect. |3.2.2| 
respectively. We consider the small mask with f^ky = 0.8 and the 
Planck-WQxPlanck-lA'i cross-spectrum simulations. Note that due 
to its low noise levels, this cross-spectrum is the most challenging 
to describe. All other cross-spectra show an even better agreement. 
The results are summarized in Eig.|Cl|and Eig.|C2| 


Figure B2. Ratio of the third order central moment {Kp for the Ax B coss 
spectrum estimator to the one from the auto-spectrum of the optimally com¬ 
bined map. The colors indicate different noise combinations according to 
the following scheme: WMAPxPlanck-10= orange, WMAPxPlanck-\00 = 
gold, WMAPxPlanck-lA3 = purple, Planck-lOxPlanck-lOO = red, Planck- 
10xPlanck-\A3 = blue, Planck-\00xPlanck-\A3 = green. The dashed line 
recalls that most of the intresting information in EE about reoinization is 
contained below f < 15. 


and where when Nf » Nf,Nf - Nf : 

Var(c/"®) - i -b NfNf) (B7) 

Var(c/) - i [2(Cff + 4NfCf + 2(Nff). (B8) 

Beyond the first term (the cosmic variance) which is dominant for 
low-f”s, the cross-spectrum picks up the noisiest of the two mea¬ 
surement while auto-spectra uses essentially the best one. 

On the other side, when both measurements have similar noise 
N'^ 

levels, Nf,N^ - —, the variances become similar: 

t ^ ^ 2 

VarCc/’"'') = ^ (2(C*)" + 2<C* -b (<)") (B9) 

Var(c/) |2(C^" -b 2<C''' -b . (BIO) 

Wether the linear term on C^dominates or not depends on the 
signal, the noise levels and the t range. As a rule-of-thumb, the 
comparision of Eq. Hzl to Eq. ( |B8^ suggests that the statistical 
loss for a cross-combination is ‘reasonable’ when the two noise 
levels are within a factor =: 3. The same kind of conclusion holds 
for the third central moment, but one can get a smaller value 
using cross-spectra for similar noise levels (Eig.|B2[l. 
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h2 1=3 



■ 0.0097 0.0729 0.1555 0.2380 - 0.0057 0.0609 0.1276 0.1942 


1=4 



■ 0.0035 0.0323 0.0680 0.1038 





Figure Cl. Normalized histograms (in black) of the 100x143 PCL estimator in the { e [2,16] range are compared to our analytic full-sky based description, 
for modell (blue) and model2 (red). The number of degree of freedom N{£) = {2£ + is reduced according to the values obtained on model 2 only 

(Fig.|4). 
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Figure C2. Normalized histograms (in black) of the 100x143 PCL estimator in the i e [2,16] range are compared to our analytical parametrization based on 
the Edgeworth expansion for model 1 (blue) and model2 (red). 
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